

% Close any parpool and open a new one
%delete(gcp('nocreate'))
%parpool(chains)


load(['Mats_T2C.mat'])
Cons_T2C.NETMISS = MISS_T2C.NETMISS;
Cons_T2C.A_MISS = MISS_T2C.A_MISS;

Chain_Out_zero = [];
if Restart == 0
    
    % Get initial estimates from multiple starting values
    Svals = 6;
    G_ij_zero = NetOut_T2C.G_ij;
    BGD_hat_s = zeros(4*Cons_T2C.dimX+3,Svals);
    A_hat_s = zeros(Cons_T2C.n,Svals);
    fval_min_s = zeros(1,Svals);
    
    parfor i=1:Svals
        i
        [NetEst_output_Svals] = Fn_EstimateNetwork_Mar2022(randn(4*Cons_T2C.dimX+3+Cons_T2C.n,1), Cons_T2C, G_ij_zero, sim_tol, 0);
        BGD_hat_s(:,i) = NetEst_output_Svals.BGD_hat;
        A_hat_s(:,i) = NetEst_output_Svals.A_hat;


    end
    fval_min_min = 1e+5
    for i=1:Svals
        if fval_min_s(1,i) < fval_min_min
            fval_min_min = fval_min_s(:,i);
            BGD_hat_min = BGD_hat_s(:,i);
            A_hat_min = A_hat_s(:,i);
        end
    end

    NetEst_output_zero = struct('BGD_hat',BGD_hat_min,'A_hat',A_hat_min,'l_M_hat',[],'beta_M_hat',[],'sigma2_M_hat',[]);
    Chain_Out_zero = struct('G_ij_b',[],'BGD_hat_b',[],'A_hat_b',[],'l_M_hat_b',[],'V_C_b',[],'beta_M_hat_b',[],'sigma2_M_hat_b',[],'fval_min_b',[],'exitflag_b',[],'toc_b',[])    

end


parfor c=1:chains
    % Start values
    if Restart==0
        Chain_Out = Chain_Out_zero;
        Chain_Out.chain_no = c;
        
        G_ij_imp = G_ij_zero;
        NetEst_output_imp = NetEst_output_zero;
        
        b_start_c = 1
        
   
    end
        

    if Restart==1
    
        [Chain_Out] = parload_b(sim_tol,c,0);
        
        b_start_c = size(Chain_Out.BGD_hat_b,2)+1;
        G_ij_imp = Chain_Out.G_ij_b(:,b_start_c-1);
        NetEst_output_imp = struct('BGD_hat',Chain_Out.BGD_hat_b(:,b_start_c-1),'A_hat',Chain_Out.A_hat_b(:,b_start_c-1),'l_M_hat',Chain_Out.l_M_hat_b(:,b_start_c-1),'beta_M_hat',Chain_Out.beta_M_hat_b(:,b_start_c-1),'sigma2_M_hat',Chain_Out.sigma2_M_hat_b(:,b_start_c-1));
        
        [c, b_start_c]
    end
    
    
    while b_start_c <= reps_chain   
        tic
        %% Estimation Step
        [NetEst_output_imp_new] = Fn_EstimateNetwork_Mar2022([NetEst_output_imp.BGD_hat; NetEst_output_imp.A_hat], Cons_T2C, G_ij_imp, sim_tol, 0);

        %% Imputation Step
        [G_ij_imp_new, A_imp_new] = Fn_ImputeNetwork_Mar2022(G_ij_imp, NetEst_output_imp_new, Cons_T2C, sim_tol, 0);


        % Reset Estimates                                
        NetEst_output_imp = NetEst_output_imp_new;
        G_ij_imp = G_ij_imp_new;

        [c, b_start_c]
        toc

        %% Save estimates and imputations
        Chain_Out.G_ij_b = [Chain_Out.G_ij_b, G_ij_imp];
        Chain_Out.BGD_hat_b = [Chain_Out.BGD_hat_b, NetEst_output_imp.BGD_hat];
        Chain_Out.A_hat_b = [Chain_Out.A_hat_b, NetEst_output_imp.A_hat];
        Chain_Out.V_C_b = [Chain_Out.V_C_b, NetEst_output_imp.V_C(:,1)];            
        Chain_Out.l_M_hat_b = [Chain_Out.l_M_hat_b, NetEst_output_imp.l_M_hat];
        Chain_Out.beta_M_hat_b = [Chain_Out.beta_M_hat_b, NetEst_output_imp.beta_M_hat];
        Chain_Out.sigma2_M_hat_b = [Chain_Out.sigma2_M_hat_b, NetEst_output_imp.sigma2_M_hat];
        Chain_Out.fval_min_b = [Chain_Out.fval_min_b, NetEst_output_imp.fval_min];
        Chain_Out.exitflag_b = [Chain_Out.exitflag_b, NetEst_output_imp.exitflag];
        Chain_Out.toc_b = [Chain_Out.toc_b, toc];

        parsave_b(Chain_Out,sim_tol,c,0);
        b_start_c = b_start_c + 1;

    end


end
